/*
 * Reaction.cpp
 *
 *  Created on: 23 May 2011
 *      Author: allan
 */

#include "Reaction.h"

// C++ includes
#include <cmath>
#include <iostream>

// FluidReact includes
#include "System.h"
#include "NumericalUtilities.h"

Reaction::Reaction()
{}

Reaction::~Reaction()
{}

void Reaction::SetNameOfReactant(string nameOfReactant)
{
	this->nameOfReactant = nameOfReactant;
}

void Reaction::SetNameOfProducts(const vector<string>& nameOfProducts)
{
	this->nameOfProducts = nameOfProducts;
}

void Reaction::SetIndexOfReactant(unsigned indexOfReactant)
{
	this->indexOfReactant = indexOfReactant;
}

void Reaction::SetIndexOfProducts(const vector<unsigned>& indexOfProducts)
{
	this->indexOfProducts = indexOfProducts;
}

void Reaction::SetStoichiometryOfProducts(const vector<double>& stoichiometryOfProducts)
{
	this->stoichiometryOfProducts = stoichiometryOfProducts;
}

void Reaction::SetEquilibriumConstant(const vector<double>& logK, const vector<double>& temperatures)
{
	this->logK = logK;
	
	this->temperatures = temperatures;
}

const unsigned Reaction::GetNumOfProducts() const
{
	return nameOfProducts.size();
}

const string Reaction::GetNameOfReactant() const
{
	return nameOfReactant;
}

const string Reaction::GetNameOfProduct(unsigned iProduct) const
{
	return nameOfProducts[iProduct];
}

const unsigned Reaction::GetIndexOfReactant() const
{
	return indexOfReactant;
}

const unsigned Reaction::GetIndexOfProduct(unsigned iProduct) const
{
	return indexOfProducts[iProduct];
}

const double Reaction::GetStoichiometryOfProduct(unsigned iProduct) const
{
	return stoichiometryOfProducts[iProduct];
}

const double Reaction::GetStoichiometryOfSpecies(unsigned iSpecies) const
{
	for(unsigned iProduct = 0; iProduct < GetNumOfProducts(); ++iProduct)
		if(GetIndexOfProduct(iProduct) == iSpecies)
			return GetStoichiometryOfProduct(iProduct);
	
	return 0.0;
}

const double Reaction::EquilibriumConstant(double T, double P) const
{
	const double log10K = LagrangeInterpolation(T, &temperatures[0], &logK[0], temperatures.size());
	
	return pow(10.0, log10K);
}






EquilibriumReaction::EquilibriumReaction()
{}

EquilibriumReaction::~EquilibriumReaction()
{}


KineticReaction::KineticReaction()
{}

KineticReaction::~KineticReaction()
{}





MineralReaction::MineralReaction()
{}

MineralReaction::~MineralReaction()
{}

void MineralReaction::AddMechanisms(const MineralReactionMechanism& mechanism)
{
	mechanisms.push_back(mechanism);
}

void MineralReaction::SetSpecificSurfaceArea(double specificSurfaceArea)
{
	this->specificSurfaceArea = specificSurfaceArea;
}

const double MineralReaction::ReactionRate(double T, double P, const System& system) const
{
	const double R = 8.314472E-3; // [R] = kJ/(mol.K)

	const double k = mechanisms[0].k;
	
	const double Ea = mechanisms[0].Ea;
	
	const double K = EquilibriumConstant(T, P);
	
	const unsigned indexOfMineral = system.GetIndexOfKineticSpecies(GetIndexOfReactant());
	
	const double A =  specificSurfaceArea * system.GetMolesOfSpecies(indexOfMineral);
	
	const double Q = system.GetKineticReactionQuotient(GetIndexOfReactant());
	
	cout << "SI = " << Q/K << "\t";
	
	return k * exp(-Ea/R * (1.0/(T + 273.15) - 1.0/298.15)) * A * (Q/K - 1.0);
}
